Polarized Wannier functions: ab-initio study of the dielectric properties of silicon and 

gallium arsenide 
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We present a first-principles calculation of the electronic properties of crystalline silicon and gal- 
lium arsenide in a uniform electric field. Polarized Wannier-Iike functions which are confined in a 
finite region are obtained by minimizing a total-energy functional which depends explicitly on the 
macroscopic polarization of the solid. The polarization charge density and the electronic dielectric 
constant are computed via finite differences. The results coincide with those of the linear response 
approach in the limit of vanishing electric field and infinite localization region. 
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The study of materials in a finite electric field is a 
difficult theoretical problem which is not definitively set- 
tled yet §. The main difficulty is that the scalar poten- 
tial which describes the field is not periodic and is not 
bounded from below: crystal momentum is no longer a 
good quantum number and no ground state exists for the 
electrons of the solid. Ab-initio studies of the response 
of crystals to uniform electric fields have been done 
mainly with linear response theory (LRT) . Dielectric con- 
stants [|[, piezoelectric tensors ||], second order non- 
linear susceptibilities [Q, and Born effective charges |^] 
of several materials have been computed with density 
functional perturbation theory (DFPT) In prac- 

tice however with these methods only the response to 
infinitesimal electric fields can be studied. 

A major issue is therefore the study of solids in a finite 
electric field. Linear as well as nonlinear susceptibili- 
ties could be extracted simultaneously from such calcu- 
lations, together with important informations on tech- 
nologically interesting phenomena, such as electromigra- 
tion, which have eluded, until now, a detailed ab-initio 
investigation Only a few direct calculations with fi- 
nite electric fields have been done with large supercells 
and artificially periodic fields 0,^. Recently Nunes and 
Vanderbilt |^ have proposed an approach to circumvent 
the difficulties associated with finite electric fields: they 
have shown that localized Wannier functions (LWFs) al- 
low one to write a functional for the energy of a solid in a 
uniform electric field. The ground state of this functional 
corresponds to a meta-stable state for the real solid , 
but, in the limit of vanishing electric field, the derivatives 
of the functional provide the same dielectric properties 
(polarization, linear j^] and non-linear Q dielectric con- 
stants) as the perturbative approach. The method has so 
far been applied only to a simple one-dimensional tight- 
binding model. In this Letter we demonstrate its appli- 
cability within a self-consistent scheme by computing the 
dielectric properties of silicon and gallium arsenide. This 



study is of relevance also for 'order N' electronic struc- 
ture calculations since LWFs are a key-ingredient for 
many of these methods. 

In a previous work [Q, we computed LWFs of sili- 
con and gallium arsenide in zero electric field within a 
self-consistent scheme in the local density approximation 
(LDA). In this work, we consider these two materials in 
a finite electric field and obtain their polarized WFs. We 
study the ability of the latter to predict the electronic 
properties of insulators in a uniform electric field and 
discuss the convergence of the dielectric properties with 
increasing size of the localization region (LR) . 

As underlined by Gonze, Ghosez, and Godby the 
validity of the Hohenberg and Kohn theorem for an 
infinite insulating crystal in a macroscopic electric field is 
questionable and in principle one should add the macro- 
scopic polarization as a basic variable. Martin and Or- 
tiz have recently reformulated the DFT for insulat- 
ing crystals to include electric field effects. However, no 
realistic polarization-dependent functional has been pro- 
posed yet and in this work we remain at the standard 
LDA level. We use equations similar to those of Ref. 0| 
except for the "exchange and correlation field", which 
stems from the polarization dependence of the exchange 
and correlation energy. 

We describe first how to implement finite electric fields 
in our scheme based on LWFs. We then show that 
the variation of the charge density induced, at first or- 
der, by a finite electric field, and computed with polarized 
LWFs, converges to the values provided by LRT. We also 
analyze the induced macroscopic polarization and show 
that, increasing the size of the LRs, its derivative with 
respect to the field converges exponentially. The asymp- 
totic values, for vanishing field, are close to the DFPT 
results. Finally we present examples of polarized WFs. 

The electronic ground state of a periodic crystal con- 
taining N interacting electrons can be described either 
by N/2 independent Bloch orbitals or by N/2 WFs. One 
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way to obtain approximate WFs for solids has been pre- 
sented in Ref. and implemented in a self-consistent 
DFT scheme in our previous work [Q. The method is 
based on the minimization of a total energy functional 
Etot [{\vo^n)Y\ which includes implicitly orthogonalization 
constraints with respect io M = Nf./2 localized orbitals 
{|'i'o,n)}, where A'^e is the number of electrons per unit 
cell, n = 1, . . . , M is the band index, and the subscript 
indicates that the orbital is centered in the unit cell 
containing the origin. The periodicity of the crystal is 
exploited by requiring that the remaining orbitals are 
obtained by translating those of the first cell by all the 
Bravais lattice vectors R;: |wi,n) = |i'o,n}- 

As already mentioned, the Hamiltonian of an infinite, 
periodic crystal in presence of a uniform static electric 
field E [|l5|, is not bounded from below. However if the 
degrees of freedom {|w^„)} are constrained to be locahzed 
in finite regions and the electric field is sufficiently small, 
the functional [01 
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(where £7 is the unit-cell volume and the macroscopic 
electric polarization P is the sum of the electronic Pe and 
ionic contributions) has a well defined minimum which 
represents a meta-stable state for the real solid [p|,p^. 
In the limit where the size of the LRs becomes infinite, 
the maximum allowed field goes to zero. For any non- 
vanishing value of the electric field E, the second term 
in Eq. ([l]) will cause a change in the periodic density 
and also a variation of the macroscopic polarization. In 
our calculations, the value of the screened electric field E 
is kept fixed, while the macroscopic polarization changes 
self-consistently. Physically we compute variations of the 
charge density, of the macroscopic polarization, and of 
the total energy with respect to the screened uniform 
electric field instead of the bare external electric field Eq. 

At zero field, the functional is minimized by almost 
orthogonal orbitals which, for sufficiently large LRs, are 
a good approximation to a set of WFs {Iuiq,™)} of the 
solid. Similar sets of WFs can also be obtained through 
a unitary transformation of the Bloch orbitals with an 
appropriate choice of the phases ||l6|. In presence of a 
finite field, the functions {Ii'o'ti)} which minimize the 
functional are still approximately orthogonal and, for 
large LRs, approximate a set of polarized WFs for the 
solid. The electronic contribution Pe to the macro- 
scopic polarization P can be related to the cen- 
ters r„ = (u'o,„|r|wo,„) of the WFs of the system by 
f2Pe = ~'^^^n=i ^n- This definition can be extended to 
non-zero values of the screened macroscopic electric field 
and to LWFs by §: 



nV, [{KE„)},E] =-2e^r„(E), 

n=l 

with r„(E) given by 
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r„(E) = }_^«Jr|<„)gt„'^.„(E), 



(3) 



; m=l 



where e is the electronic charge, and (5j;f„(E) — 
2Si^oSrn,n ^ (^^Pmbo'n) thc first-ordcr approximation to 
the inverse of the overlap matrix between LWFs. 

This extended definition of Pg allows one to study the 
dielectric properties of insulating crystals in finite electric 
fields. The validity of this definition rests on the agree- 
ment of the values of the dielectric properties computed 
in the limit of vanishing field with the LRT results . 
In this limit the equivalence between the two methods 
can be demonstrated analytically. The electronic dielec- 
tric tensor (E) is related to the macroscopic electronic 
polarization Pg by 
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where a and (3 indicate Cartesian coordinates. It can be 
obtained by finite differences of thc macroscopic polar- 
ization Pe(E) for different values of the screened electric 
field E. Pe(E) refers to Pe [{|dJ;„)},E] of Eq. (|) com- 
puted for the orbitals {|w^„)} which correspond, for a 
given electric field E, to the minimum of thc functional 
F . 

Thc search for the minimum of the total-energy func- 
tional F requires its functional derivative with respect 
to thc variational degrees of freedom (r|ti^„). Differen- 
tiating Eq. (|l]) with respect to (r|w^„), the polarization 
term gives an additional contribution with respect to the 
zero-field case, which can be recast as: 
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The LWFs have been represented on a uniform real- 
space mesh and (rlwo'n) is allowed to be non-zero only 
inside cubic regions with size 2a lr. As in the zero-field 
case, the localization of the LWFs makes the summa- 
tions in Eqs. (|^) and (||) finite. The periodic valence 
charge density n(r; E) (whose expression is the same 
as in the zero-field case) is a sum of localized contri- 
butions. The Fourier components n(G; E) allow us to 
compute the Coulomb electrostatic energy per unit cell 
as Eh = §EG^o^''i*(G;E)Vff(G;E), where Vh(G;E) 
is the Fourier transform of the Hartree potential. The 
G = component of the Hartree potential is included 
in the screened electric field E and it is not explicitly 
evaluated. 

We have applied this approach to crystalline silicon 
and gallium arsenide. The unit cell contains 2 atoms and 
the electronic ground-state can be described by M = 4 
LWFs which are confined within LRs centered on the four 
bonds in the unit cell at the origin. Our calculations have 
been carried out within DFT-LDA and with norm- 
conserving pseudopotentials in the Kleinman-Bylander 
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form All technical details and parameters are the 

same as those given in Ref. . 

In Fig. 1^, we show the charge density induced by a 
field in the [100] direction and with an intensity equal 
to e |E| = 10-3 a.u. « 5.14 x 10^ V/m. The induced 
valence charge density An(r) — n{r; E) — n(r, 0) com- 
puted with LWFs obtained for different values of qlr 
is compared in Fig. [l](a) with the result of the conven- 
tional linear response approach [^ n^^)(r) = E ^^g^ 

E— 

obtained using a Bloch orbital representation of the elec- 
tronic wave-functions ||2^ . The weak intensity of the field 
ensures a linear behaviour for both silicon and gallium ar- 
senide. Fig. |^(a) shows that all the features of the charge 
density linearly induced by the field are correctly repro- 
duced with small LRs. Localization regions containing 
216 atoms {uLR/a = 34/24) are large enough to obtain 
an induced charge density which is practically undistin- 
guishable from the LRT result on the scale of the figure. 
The difference is everywhere less than 4.7 x IQ-^ e/cell for 
silicon and 6.1 x IQ-^ e/cell for gallium arsenide. In the 
case of gallium arsenide, this error is smaller than that 
associated with the use of a semi-local pseudopotential 
instead of the Kleinman-Bylander form used here. 

In Fig. H, we display the values of the high-frequency- 
dielectric constants of silicon and gallium arsenide com- 
puted with our method at E = for different sizes of 
the LRs. Fitting the data with an exponential function 
e°° — Aexp{—{aLji/a) /a) + B, we obtain in the limit 
aLR ^ oo the values B — 13.4 ± 0.2 for silicon and 
B = 11.6 ± 0.1 for gallium arsenide. The LRT results are 
12.9 and 11.4 respectively, thus, supporting the use of 
LWFs in the study of the dielectric properties of sohds. 
The convergence of the dielectric constant with increas- 
ing values of qlr^ is slower than that of structural prop- 
erties The exponential fitting curves correspond to 
the parameters: A = -13.2 ± 0.3, a = 0.81 ± 0.05 for 
silicon and A = -10.4 ± 0.2, a = 0.83 ± 0.05 for gal- 
lium arsenide. For the largest localization size consid- 
ered {(ilr/o. = 2 which contains 512 atoms) we obtain 
e°° = 12.3 for silicon and e°° — 10.6 for gallium arsenide, 
which underestimate the LRT results by 5% and 7% re- 
spectively. 

In Fig. |3|(a), we show the polarized LWF computed 
at the theoretical lattice constant, a^R/a — 34/24, and 
with an electric field e |E| — 10^^ Ha/a.u. in the [100] 
direction for silicon and gallium arsenide. The electric 
field displaces the center of each LWF with respect to 
the zero- field case, as indicated with arrows in Fig. §(b). 
By summing over the four functions within the unit cell 
at the origin, the individual displacements perpendicu- 
lar to the field cancel out and the induced macroscopic 
polarization is parallel to the field, as expected for solids 
with cubic symmetry. 

In conclusion we have shown that LWFs, a key ingredi- 
ent of several 'order N' methods for electronic structure 



calculations, also allow one to study solids in presence of 
a finite electric field. The induced charge density and the 
dielectric constant obtained from finite differences con- 
verge, in the limit of large LRs and small electric fields, 
towards the LRT values. Furthermore a first example of 
approximate, polarized WFs has been provided. 
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FIG. 1. Valence charge density An(r) (in units of e/cell) 
induced by a macroscopic electric field E in the [100] direction 
(with e |E| = 10""^ a.u.). Panel (a): Comparison of results 
obtained for two different values of fliij/a (the number of 
atoms inside a LR is in square brackets) with the LRT result. 
The induced density is represented along the [111] direction 
at the theoretical unit cell volume and d = \/^a. Panel (b): 
Contour-plots in the (Oil) plane of An(r) for the largest LR 
size. The contours are drawn at constant intervals of 0.02 
e/cell. The thicker contour denotes An(r) = 0, the continuous 
(dashed) contours are for An(r) > (An(r) < 0). 



FIG. 2. Linear dielectric constant e°° of silicon and gal- 
lium arsenide computed with LWFs using different LR sizes. 
The values of e°° for ulr —> oo axe obtained by numerical 
extrapolation, e^m- indicates the LRT value. Calculations 
have been performed at the theoretical lattice constants 
ao = 10.20 a.u. for silicon and ao = 10.48 a.u. for gallium 
arsenide. 

FIG. 3. Panel (a): Contour-plots in the (Oil) plane of 
the polarized LWFs in a.u. obtained for an electric field 
e |E| — 10~^ a.u. oriented along [100]. These polarized 
LWFs are constrained to be zero outside a LR with size 
dLR/a, = 34/24 and centered on the bond oriented along the 
[111] direction. The contour intervals are 0.01 a.u. for neg- 
ative values (dashed lines) and 0.1 a.u. for positive values 
(continuous lines). Panel (b): Contour-plots of the variation 
of the LWFs (r|At;o.„^i) = {r|u^„=i) - {r\vf^^Zi). The con- 
tour interval is 10"'' a.u. for silicon and 2 x 10"'* a.u. for 
gallium arsenide. The LWF centers are indicated by empty 
circles while their variation is depicted by arrows whose length 
has been magnified by a factor 100. 
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